"""
Common utilities and base functions for differential operators across all coordinate systems.
"""

using LinearAlgebra
using SparseArrays
using FFTW

# Export main operator creation functions
export OperatorType, BoundaryCondition, make_general_poisson_solver
export uniform_discretization, get_operator_stencil
export apply_boundary_conditions!, cache_operator

"""
Type alias for operator functions
"""
const OperatorType = Function

"""
Abstract base type for boundary conditions
"""
abstract type BoundaryCondition end

"""
Dirichlet boundary condition: u = value
"""
struct DirichletBC <: BoundaryCondition
    value::Float64
end

"""
Neumann boundary condition: ∂u/∂n = value
"""
struct NeumannBC <: BoundaryCondition
    value::Float64
end

"""
Robin boundary condition: au + b∂u/∂n = value
"""
struct RobinBC <: BoundaryCondition
    a::Float64
    b::Float64
    value::Float64
end

"""
Periodic boundary condition
"""
struct PeriodicBC <: BoundaryCondition
end

"""
    uniform_discretization(grid) -> Vector{Float64}

Get uniform discretization spacing for each dimension of the grid.
"""
function uniform_discretization(grid)
    return [spacing(grid)[i] for i in 1:ndim(grid)]
end

"""
    get_operator_stencil(order::Int, dim::Int) -> Tuple{Vector{Int}, Vector{Float64}}

Get finite difference stencil coefficients for given order and dimension.
"""
function get_operator_stencil(order::Int)
    if order == 1
        # First derivative central difference
        return ([-1, 1], [-0.5, 0.5])
    elseif order == 2
        # Second derivative central difference  
        return ([-1, 0, 1], [1.0, -2.0, 1.0])
    elseif order == 4
        # Fourth-order second derivative
        return ([-2, -1, 0, 1, 2], [-1/12, 4/3, -5/2, 4/3, -1/12])
    else
        error("Unsupported stencil order: $order")
    end
end

"""
    apply_boundary_conditions!(matrix, vector, bcs, grid)

Apply boundary conditions to operator matrix and vector.
"""
function apply_boundary_conditions!(matrix, vector, bcs, grid)
    # This is a simplified implementation - full implementation would
    # handle different BC types and grid geometries
    dims = size(grid)
    total_points = prod(dims)
    
    for (dim_idx, bc) in enumerate(bcs)
        if isa(bc, DirichletBC)
            # Apply Dirichlet BC: set diagonal to 1, off-diagonals to 0
            # This is a simplified version
            continue
        elseif isa(bc, NeumannBC)
            # Apply Neumann BC: modify stencil at boundaries
            continue
        elseif isa(bc, PeriodicBC)
            # Connect boundary points
            continue
        end
    end
    
    return matrix, vector
end

"""
    cache_operator(grid, operator_func, cache_key::Symbol)

Cache an operator for reuse to improve performance.
"""
function cache_operator(grid, operator_func, cache_key::Symbol)
    # Simple caching mechanism - could be enhanced with LRU cache
    if !hasfield(typeof(grid), :_operator_cache)
        # Add cache field if not present
        grid._operator_cache = Dict{Symbol, Any}()
    end
    
    if haskey(grid._operator_cache, cache_key)
        return grid._operator_cache[cache_key]
    else
        operator = operator_func()
        grid._operator_cache[cache_key] = operator
        return operator
    end
end

"""
    make_general_poisson_solver(grid, bcs) -> Function

Create a general Poisson equation solver for any grid type.
"""
function make_general_poisson_solver(grid, bcs)
    L_matrix, L_vector = make_laplace_matrix(grid, bcs)
    
    # Use sparse solver
    function poisson_solver(rhs, out)
        # Solve L*u = rhs
        solution = L_matrix \ (vec(rhs) - L_vector)
        out .= reshape(solution, size(rhs))
        return out
    end
    
    return poisson_solver
end

"""
    make_laplace_matrix(grid, bcs) -> Tuple{SparseMatrixCSC, Vector}

Create Laplacian matrix and boundary condition vector for given grid and boundary conditions.
This is a dispatch function that calls grid-specific implementations.
"""
function make_laplace_matrix(grid, bcs)
    error("make_laplace_matrix not implemented for grid type $(typeof(grid))")
end

"""
    validate_operator_accuracy(operator_func, analytical_func, grid, tolerance=1e-10)

Validate numerical operator accuracy against analytical solution.
"""
function validate_operator_accuracy(operator_func, analytical_func, grid, tolerance=1e-10)
    # Create test function values on grid
    test_values = zeros(size(grid))
    for (i, coord) in enumerate(grid_coordinates(grid))
        test_values[i] = analytical_func(coord...)
    end
    
    # Apply numerical operator
    numerical_result = operator_func(test_values)
    
    # Compare with analytical result
    analytical_result = zeros(size(grid))
    for (i, coord) in enumerate(grid_coordinates(grid))
        analytical_result[i] = analytical_func(coord...)  # Should be analytical derivative
    end
    
    error_norm = norm(numerical_result - analytical_result, Inf)
    return error_norm < tolerance
end

"""
    grid_coordinates(grid)

Get iterator over all coordinate points in the grid.
"""
function grid_coordinates(grid)
    # This would need to be implemented for each grid type
    error("grid_coordinates not implemented for grid type $(typeof(grid))")
end

"""
    conservative_stencil_weights(r_inner, r_outer, volume)

Calculate conservative finite difference weights for spherical/cylindrical coordinates.
"""
function conservative_stencil_weights(r_inner, r_outer, volume)
    # For conservative operators in curved coordinates
    factor_inner = r_inner^2 / volume
    factor_outer = r_outer^2 / volume
    return factor_inner, factor_outer
end

"""
    periodic_extension!(array, grid)

Extend array with periodic boundary conditions for FFT-based operators.
"""
function periodic_extension!(array, grid)
    # Implement periodic extension for FFT operators
    return array
end

"""
    fft_frequencies(n, L)

Generate FFT frequency array for spectral operators.
"""
function fft_frequencies(n, L)
    # Generate k-space frequencies for FFT-based operators
    k = fftfreq(n, 2π/L)
    return k
end